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ABSTRACT 

We study the formation of a supermassive black hole (SMBH) binary and the shrinking 
of the separation of the two holes to sub-pc scales starting from a realistic major merger 
between two gas-rich spiral galaxies with mass comparable to our Milky Way. The 
simulations are the first of this kind carried out with an Adaptive Mesh Refinement 
(AMR) code, in this case the RAMSES code, and the first capable to resolve separations 
as small as 0.1 pc. The collision of the two galaxies produces a gravo-turbulent rotating 
nuclear disk with mass (~ 10^ M©) and size (~ 60pc) in excellent agreeme nt with 
previous S PH simulations with particle splitting that used a similar setup (|Mayer| 
et al.|2QQ7| but were limited to separations of a few parsecs. The AMR results confirm 
that the two black holes sink rapidly as a result of dynamical friction onto the gaseous 
background, reaching a separation of 1 pc in less than 10^ yr. We show that the 
dynamical friction wake is well resolved by our model and we find good agreement 
with analytical predictions of the drag force as a function of the Mach number. Below 
1 pc, black hole pairing slows down significantly, as the relative velocity between the 
sinking SMBH becomes highly subsonic and the mass contained within their orbit 
falls below the mass of the binary itself, rendering dynamical friction ineffective. In 
this final stage, the black holes have not opened a gap as the gaseous background is 
highly pressurized in the center. Non-axisymmetric gas torques do not arise to restart 
sinking in absence of efficient dynamical friction, at variance with previous calculations 
using idealized equilibrium nuclear disk models. We believe that the rather "hot" 
Equation-of-State we used to model the multiphase turbulent ISM in the nuclear 
region is playing an important role in preventing efficient SMBH sinking inside the 
central parsec. We conclude with a discussion of the way forward to address sinking 
in gaseous backgrounds at sub-pc scales approaching the gravitational wave regime. 
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1 INTRODUCTION 

The coalescence of two supermassive black holes (SMBHs), 
namely black holes in the mass range 10^ - 10^ Mo, would 
produce the loudest gravitational wave signals in the Uni- 
verse, and is the main target of planned low-frequency grav- 
itational wave experiments, such as space-based laser inter- 



ferometers and pulsar timing arrays JVecchio 2004; Sesana 



[et al.||2009| . Since all massive galaxies are believed to host 
a central supermassive black hole, in the current cosmolog- 
ical paradigm in which galaxies grow via repeated mergers, 
mergers between such black holes could be common, espe- 
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cially at high redshift when the galaxy merger rate is higher. 
However, the merger rate of SMBHs does not follow trivially 
from that of their host galaxies. Once the two galaxy cores 
have merged, leaving no distinct substructure at hundred of 
parsecs/kiloparsec scales, the two black holes will need to 
reduce their separation to less than 0.01 pc before they can 
start to lose orbital energy efficiently via gravitational wave 
emission and eventually coalesce ( ,Begelman, Blandford fc] 
Rees||198"0 ). Loss of orbital energy can occur via dynamical 



friction onto the stellar background (Milosavljevic & Merritt 
|2001|) or due to the gas drag { Esc ala et al.|2004 ) . Both mech- 
anisms are relevant since SMBHs are inferred to exist at the 
center of both gas-rich spirals and gas-poor ellipticals/SOs 
( jVolonteri, Haardt k Gultekin|2008t . The first one is known 
to become ineffective when the binary begins to harden; at 
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this stage the 3-body interactions between the binary and 
individual stars deplete its loss cone, which can be refilled 
only via some orbital diffusion mechanism that brings fresh 
stellar material from other regions of the galaxy (Berczik^ 
eraL][2005t . 



In recent years it has been shown that SMBH decay 
in a gas dominated system is much more effective at form- 
ing rapidly a binary of SMBHs, in only a few million years 
following the major merger of two moderately gas-rich disk 
galaxies ( [Mayer et al.|2007| payer Kazantzidis|2008t . The 
reason is that the decay occurs in a much denser medium rel- 
ative to the stellar dominated case, and because the black 
holes move with Mach numbers of order unity or slightly 
larger, the dynamical friction onto a gaseous medium is ex- 
pected to be more effective than in a collisionless system 
(Ostriker 1999). Indeed the two black holes spiral down 



to parsec scales in a gaseous, dense nuclear disk a hun- 
dred pc in size formed by the dramatic gas inflow in the 
merger. Such nuclear disks have been found in high resolu- 
tion multi-wavelength observations of nearby merger rem- 
nants ( Downes Solomon|1998||Davies et al.|200 7', etc.). In 
minor mergers efficient black hole pairing at scales of tens 
of parsecs, namely even before a bound binary can form, re- 



quires high gas fractions (> 30%) in both galaxies (Callegari 
et al. 12009) 2011| strengthening even further the crucial role 



of gas in the pairing and merging of SMBHs. Other inves- 
tigations starting when the two black holes form already a 
loosely bound pair in a nuclear disk have suggested that the 
decay can continue to sub-pc scales under certain conditions 
but have utilized pre-defined models of nuclear disks (Escala| 
[et al.|2005||Dotti et al . 2006, 2007; Quadra et aT7|2009| ) rather 
than starting from a realistic galaxy merger and appeal to 
tidal torques rather than dynamical friction when the binary 
begins to harden. At the smallest scales they have assumed 
very specific configurations such as that of two black holes 
already at less than 0.1 pc separation embedded in a gap 
within the disk ( Cuadra et al.|[2009 ). Gap opening implies 
the transition to a regime in which tidal torques from the 
surrounding disk become the dominant mode to extract en- 
ergy and angular momentum from the binary. However, the 
configuration of the host on which the two black holes are 
found at less than parsec scales is not really known because 
no computation exists that can reach such a stage starting 
from a realistic merger, and consequently the decay in such 
regime is not yet explored by three-dimensional simulations. 

Until now, studies of binary black hole formation and 
shrinking in a gaseous environment have been based exclu- 
sively on SPH simulations, which might be unable to cap- 
ture the turbulent nature of the fiow in the disk. Turbulence 
arises as a result of shocks during the final galaxy collision, 
and later in the nuclear disk as a result of its gravitational 
instability ( ^Mayer et al.||2007| [Mayer &: Kazantzidis||2008 ) . 
Furthermore, dynamical friction, which is the central phys- 
ical process involved, is numerically challenging to capture 
since it involves both the effect of the local overdensity trail- 
ing the black holes and that of the larger scale torques and 
tidal distortions generated by the surrounding mass distri- 
bution ( ,Colpi et al. J999 ) 



While TreeSPH codes such as GASOLINE (Wadsley et al 



2004) are typically well suited to address processes in the 



resolve gradients between media of different densities (such 
as the overdense wake and the surrounding background), 
may cast doubts on the quantitative results regarding the 
strength of dynamical friction. More generally, confirmation 
of the results of SPH calculations on the effectiveness of the 
SMBH binary formation process in a gaseous environment 
is desirable. Due to the many scales involved, adaptive mesh 
refinement (AMR) simulations are ideally suited to the prob- 
lem. In this paper, we used the RAMSES code ( Teyssier|2002 ) 
to follow both the galaxy merger at large scale and the bi- 
nary SMBHs merger at small separation. RAMSES is based on 
a Particle-Mesh N-body integrator for dark matter and stars 
and a second-order unsplit Godunov scheme for the gas. Its 
shock capturing abilities and low intrinsic numerical viscos- 
ity (we used the HLLC Riemann solver and the MinMod 
slope limiter) make it ideal to address the problem. It also 
features a multigrid Poisson solver that is both fast (of or- 
der N) and second-order accurate ( Guillet Teyssier [2011 ), 
this being another desirable feature for the problem at hand. 
Finally, the super-Lagrangian refinement strategy we have 
used in this project allows in principle to reach smaller sepa- 
rations than possible with an SPH code in a realistic galaxy 
merger at a reasonable computational cost. 

In this paper we present the first AMR simulations of 
the formation and hardening of a binary of SMBHs. These 
are also the first simulations of this kind that, starting from 
a realistic galaxy merger, can follow the decay of the two 
SMBHs to a separation as small as 0.1 pc. 

The paper is organized as follows. In section [2] we de- 
scribe the hydrodynamical simulations we performed. Re- 
sults on SMBHs pairing and dynamical friction due to gas 
are presented in section [Sj Section [4] is devoted to discussion 
and conclusions. 



2 HYDRODYNAMICAL SIMULATION 
PARAMETERS 

We use the AMR hydrodynamical code RAMSES fTeyssierl 



2002) to model the evolution of a galactic major merger in 



which each galaxy hosts a SMBH at its center. The galac- 
tic model, inspired by the one from [Mayer Kazantzidis| 
(2008), consists in a NFW isotropic dark matter halo, a ro- 
tating exponential stellar disc, a Hernquist bulge, a thin ro- 
tating exponential gas disc and a single motionless and non- 
accreting SMBH particle placed at the center of the galaxy. 

The dark matter halo has a virial rotation velocity 
'^^200 = 138 km/s (M200 = 8.7 x lO"*^"*^ M©), a concentration 
parameter c = 9.0, a dimensionless spin parameter A = 0.05 
and is made of Ndm 7 x 10^ particles. The disk has a mass 
Md = 0.05 M200, a scale length Rd = 3.6 kpc, a thickness 
Hd = 0.1 Rd and is made of A^* = 1.8 x 10^ particles. The 
bulge has a mass Mb = 0.02 M200, a scale radius a = 0.1 Rd 
and is made of Nb = 8x 10^ particles. As in |Kazantzidis et al. 



(2005), the mass of the SMBH particle Mbh = 2.5 x 10^ M(: 



is chosen consistently with the Mbh — cr relation ( Ferrarese 
k Merritt||2000| [Tremaine et al.| f2002). The gaseous disk is 



domain of self- gravitating systems, to which category dy- 
namical friction belongs, particle noise and the difficulty to 



initialized on the AMR grid as a continuous density field 
with a total mass Mg 4.18 X 10^ M©. 

For the merger model, we use the same parabolic copla- 
nar orbit presented in Mayer & Kazantzidis (2008). During 



the interaction, both SMBH particles keep at the center the 
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galactic cores in which they were embedded. After a few 
pericentric passages, dynamical friction finally makes the 
two galaxies merge at t :^ 5.244 Gyr. The first 4 Gyr of 
the simulation have been run at low resolution to let us fo- 
cus on the final merger of the galactic cores at much higher 
resolution. 

To avoid binary relaxation among DM/stellar particles, 
the gravitational softening of every particle except for the 
two SMBH particles cannot be smaller than tmin = 3 pc, 
consistently with our particle mass resolution. Meanwhile, 
the two SMBH particle dynamics and the gas dynamics fol- 
low the local resolution of the AMR grid without any limi- 
tation. 

2.1 Gas physics and equation-of-state 

We adopt in this paper the following thermodynamical 
model. It is based on a polytropic equation of state sim- 



ilar to [Teyssier et al. (2010) and Bournaud et al. (2010). 
For the hot virialized halo, defined as tih < Tih with Uh = 
10"^ H/cc), we use a polytrope with T = (uH/rih)^^'^ 
with Th = 4 X 10^ K and Th = 5/3. For the disk, defined 
as uh > Tih, we use first an isothermal equation-of-state 
with T = 10^ K. At the end of the merger, the two black 
holes are embedded in the nuclear gas core of the remnant 
galaxy. The thermodynamics of this inner region is the key 
parameter of such a study. It is interesting to investigate the 
impact of various models on the black hole orbital decay. 
First, we assumed that a strong heating source such as an 
active galactic nucleus could completely shuts down radia- 
tive cooling in the high-density regions, resulting in a purely 
adiabatic gas (Fa = 5/3). Spaans Silk| ( [200 ) showed that 
the thermodynamic state of a solar metallicity gas heated by 
a starburst can be well approximated by an ideal gas with 
adiabatic index Fa = 1.3 — 1.4 over a wide range of densities. 
Consequently, we also adopted a second model for which the 
high-density gas is more dissipative than in the first model 
and follow a Fa = 7/5 polytropic equation-of-state. For both 
models, using as characteristic density Ud = 1000 H/cc, we 
use the following equation-of-state across the whole disk 



[lO" K,Td{nH/ndf^-'] 



(1) 



As shown by [Escala et al] ( |2005t and |Dotti et al] ( |2007| ), 

a sound speed in the range 60 — 80 km/s, corresponds to a 
pressure scale height comparable to that deduced from ob- 
servations of circumnuclear disks in nearby galaxies ( [Downes' 
Solomon] 19981 . In our model, we adopted a characteristic 
temperature Td = 3 x 10^ K or equivalently a characteristic 
sound speed Cs,d = 75 km/s. 



3 RESULTS 

A few Myr before the final merger of the two gaseous 
cores, during the last pericenter, the symmetry of the gas 
distribution changes dramatically. In Figure [l] one can see 
that the gas distribution, perfectly symmetric before the last 
pericenter, becomes asymmetric afterwards. This is due to 
small scale turbulence between the two galactic cores (top 
left view) that made the shocks within the gas asymmet- 
ric during the last pericenter. Afterwards, one of the galac- 
tic core (and SMBH within) ends up with a bit slower and 
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]= 5.2443 Gyr 5.2500 Gyr 



Figure 1. Mass- weighted gas density maps during last pericenter 
and final merger. The line-of-sight is perpendicular to the orbital 
plane and the maps are 1.8 kpc wide. While fairly symmetric 
before the pericenter (top left), the density distribution becomes 
clearly asymmetric after the pericenter (middle). After the final 
merger, a thick ~ 10^ Mq gaseous nuclear disk is formed (bot- 
tom), in which the two SMBHs start orbiting. 



with a little more mass than the other (bottom right view). 
An asymmetrical SMBH injection in the nuclear disk re- 
sults from this symmetry breaking : one of the black hole 
is injected a few parsecs away from the nuclear disk center, 
while the other starts orbiting with an apocenter as large as 
r - 100 pc. 



3.1 Nuclear disk formation and turbulence 
dissipation 

After an additional 4 Myr , the two galactic cores finally 
merge together. In the Fa = 5/3 case, the core coalescence 
leads to the formation of a 140 pc thick gas spheroid while 
in the Fa = 7/5 case, the nuclear region is more disky with 
a disk thickness of ^ 60 pc. In both cases, the enclosed gas 
mass within 100 pc is ^ 10^ M©. The parameters of the 
nuclear disk we obtain from this galaxy merger simulation 
are consistent with nuclear regions observed in Ultra Lumi- 
nous Infrared Galaxies (ULIRGs; Downes &: SolomQii||1998 



Sanders Mirabel||1996 ), except that our model is missing 
the intense star formation processes of the observed nuclear 




Figure 2. 200 pc-wide nuclear disk (Fa = 7/5 case) at t(Gyr) = 5.25 (left), 5.256 (middle) 5.2571 (right). From top to bottom are shown 
the Toomre parameter, velocity dispersion, sound speed, face-on and edge-on gas surface density maps. 
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regions. Figure [2] shows the velocity dispersion, sound speed, 
surface density and Toomre stabihty parameter maps of this 
disk (Fd = 7/5 case). The Toomre parameter is above 1 
across the whole disk. The velocity dispersion cTv is domi- 
nated by the sound speed c^. 

Compared to SPH calculations made by [Mayer Ez\ 
Kazantzidis ( 2008 ) with the same merger model, the overall 



properties of the disk are qualitatively in good agreement 
with our AMR simulations. The nuclear disk formed in the 
SPH simulation was however found to be marginally grav- 
itationally unstable. This caused the formation of a strong 
spiral wave and the collapse of the disk as the angular mo- 
mentum was driven outward. The resulting inflow was at the 
origin of a fast decrease of the ambient gas density around 
the SMBH and of the stalling of the orbital decay. Testing 
the robustness of this result using a different thermodynam- 
ical model and a different numerical method is precisely the 
main justification of the present work. A more stable nuclear 
disk configuration would result in a larger gas density and 
in a stronger drag force. In our AMR model, without any 
other source of turbulence than gravity, the gas velocity dis- 
persion is quite high at the beginning, and a rather strong 
spiral mode can be seen in the density maps. But both the 
spiral wave and the turbulence are slowly dissipated from 
(jv — 70 km/s to (Jv — 30 km/s (Fig. |2] second row) in 
the nuclear disk, over a timescale tcross — h/ av — 2 Myr. In 
contrast to the SPH results, we do not see a sustained trans- 
port of angular momentum outward, and the disk settles in 
a stable pressure equilibrium, with high gas density but also 
high gas sound speed (cs ~ 200 km/s). 
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Figure 3. SMBH particles relative separation evolution for the 
low-resolution runs (Axmin = 3 pc) with Fa = 5/3 (red) and 
Fd = 7/5 (blue) and the high-resolution run (Z\X 777,272 — 0.1 pc) 
with Fd = 7/5 (black). The blue (resp. black) dashed line corre- 
sponds to the low (resp. high) spatial resolution limit. The red 
dashed line corresponds to the black hole separation limit where 



the binary becomes hard (Mb 



: 2Mb 



: M(r < 3 pc)). 



3.2 SMBHs orbital decay 

We follow the evolution of the two SMBH particles in this 
nuclear region for the different models. Figure |3] shows the 
evolution of the black hole relative separation. The low reso- 
lution runs are in good agreement with [Mayer Kazantzidis] 
(20081 results. While the hot (Fd = 5/3) thermodynamical 



model leads to a stalling of the black hole orbital decay 
around r ~ 40 pc, the cold model (Fd = 7/5) let the black 
hole relative separation falls down to the numerical resolu- 
tion (Ax = 3 pc) after only 30 Myr. These first experiments 
are however not conclusive, since, as we will demonstrate in 
the next section, this resolution is too low to resolve properly 
the wake causing the hydrodynamical friction in the nuclear 
disk. 

In our high resolution run, shown as the black line in 
Figure [3] the environment of the SMBH is much better re- 
solved. As a consequence, the SMBH relative separation de- 
crease from 100 pc down to 1 pc in less than 10 Myr. Al- 
though the black hole binary system reach very briefly a rel- 
ative separation close to our resolution limit (Ax = 0.1 pc), 
the orbital separation of the black holes is stalling well above 
the resolution limit and settles around 2 parsecs where the 
binary system becomes hard (2Mbh = M(r < 3 pc), red 
dashed line). Using a different thermodynamical model than 
the previous SPH simulations results in a denser, more sta- 
ble nuclear disk, but, as we now show in more details, this 
also leads to inefficient hydrodynamical friction and failure 
of the model to harden the binary system down to sub-parsec 
scales in the center of this nuclear disk. Note however that, 
in our case, the dynamical friction time scale is increased 



by roughly one order of magnitude, so that orbital decay is 
not stopped entirely and SMBH pairing will probably take 
place after several tens of Myr. We couldn't follow the late 
time evolution over such a long period by lack of sufficient 
computational resources. 



3.3 Dynamical friction in a gaseous medium 

Figure |4] shows the gas overdensity induced by a SMBH 
particle during an orbit where the black hole is about r = 
20 pc from the center of the nuclear disk. On each panel, the 
position of the black hole is marked by a black dot and the 
orientation is taken so that the relative velocity of the black 
hole in its surrounding gas is pointing rightward. The black 
hole induces a shock and a trailing hydrodynamical wake, 
which both get stronger as the black hole reach a transonic 
regime where the Mach number is defined using the relative 
velocity of the SMBH with respect to the gas 



M 



s(Xbh) 



Cs(Xbh) 



(2) 



The hydrodynamical wake exerts a gravitational drag on 
the black hole proportional to 47Tpc^R^YL^ w here the Bond i 
radius is expressed as Rbh = GMbu/c^ (Ruffert 1996). 



The true efficiency of the dynamical friction exerted by the 
gaseous medium on the black hole can be expressed by the 
dimensionless correction factor 



(gas) 



47rp(GMBH/cs)2 



(3) 
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Figure 4. Gas Overdensity maps behind the SMBH particle (the relative velocity of the SMBH compared to the surrounding gas is 
pointing rightward). The hydrodynamical wakes and Mach cones are shown for increasing values of A4 = Vbh/gas/cs- These wakes are 
clearly stronger and make the dynamical friction much more efficient when the black hole is in a transonic regime (A4 = 1.14, 1.18, 1.23). 



The analytical study by Ostriker ( 19991 provided a formula 

for /^^"^) 



/>(gas) 
J subsonic 



/.(gas) _ 
J supersonic 



1 

1 



1, fl + M 



M 



In {M^ 



1) +lnA 



(4) 
(5) 



where the Coulomb logarithm In A = ln(rmax/Tmin) accounts 
for the maximum and minimum radial contributions to the 
drag, as in the standard case of dynamical friction exerted 
by a stellar background ( Chandrasekhar 1943). Note that 



previous numerical studies have claimed that these analyti- 
cal formulae were probably overestimating the effect of dy- 
namical friction for the transonic regime (Escala et al.|2004 



Sanchez- Salcedo Brandenburg|2001 ). In order to compare 



our numerical drag to these analytical and numerical esti- 
mates, we computed the gravitational force of the perturbed 
density field seen by the orbiting black hole particle. FigurejS] 
shows the dynamical friction dimensionless factor f^^'^''\ The 
numerical values obtained during one orbit of the distant 
black hole particle is in good agreement with the analytical 
prediction of |Ostriker| (|1 999| if we use In A 3. Deviations 
from the analytical model could be easily explained by the 
time-dependent nature of the SMBH orbit, while the the- 
ory of Ostrikerl (11999 ) is based on a strictly stationary flow 



around the black hole. 

In Figure [6] we plot the fractional contribution to the 
total gravitational drag of plan-parallel slabs of gas perpen- 
dicular to the propagation axis of the BH. In this plot, the 
x-coordinate is the distance of each slab (of size 3 pc) to the 
BH. We see immediately that for various Mach numbers, 
the maximum radius that contributes to the drag force is 
roughly rmax — 2.5 pc. Moreover, the fractional drag pro- 
file is well resolved by the cell size of our simulation, which 
quite naturally corresponds to the minimum scale that con- 
tributes to the drag rmin — 0.1 pc. From these two numbers. 
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Figure 5. Dimensionless factor /(s^^) of the dynamical friction 
force as a function of the Mach number M. = '^^b/t/gas/^s- An- 
alytical prediction by [Ostriker] ( |1999 ) for a Coulomb logarithm 
InA = 3 (red) and InA = 3.5 (green) compared to numerical re- 
sults from [ Escala et ar] ( |2004| ) (blue) and this work (black dots), 
corresponding to a numerical Coulomb logarithm InA ~ 3.2. 



we can estimate the Coulomb logarithm in our simulation as 
InA = Inrmax/Tmin — lu 25 ^ 3.2. We see in Figure [s] that 



the analytical model of [Ostriker (1999) using InA = 3 (red) 



or In A = 3.5 (green) is a perfectly good fit to our numerical 
data. From Figure |6] we see that resolving the SMBH envi- 
ronment with sub-parsec resolution is mandatory in order to 
resolve the wake properly. At Mach number slightly above 
unity, the wake profile sharply declines towards the BH po- 
sition, and even higher resolutions would probably result in 
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Figure 6. Spatial contribution to the gravitational drag the gas 
exerts on the black hole particle at different Mach number. The 
strength of the drag peaks around transonic regime {M. = 1.13, 
red). Whatever the Mach number, most of the gravitational drag 
comes from the 2 — 3 pc region trailing behind the black hole. 



a stronger drag force, by effectively raising the value of the 
Coulomb logarithm. 

3.4 A transition from fast to slow orbital decay 

Let us consider a simple model of orbital decay. Assuming 
that the black hole evolves on circular orbits with circular 
velocity Vc\vc{r) and undergoes a hydrodynamical friction 
force opposite to its velocity vector and of amplitude Fdf- 
The BH angular momentum writes 



L — Munrvcirc 

and the angular momentum losses write 



|r X Fd 



We can derive the orbital decay timescale as 
L ~ i^df 



Tod 



(6) 



(7) 



(8) 



When the BH is evolving at large radii, say 20 pc, the gas 
density is rather low, p = 6 x 10^ H/cc, as is the gas sound 
speed with Cs — 230 km/s. The Mach number of the BH 
relative speed is however measured in the simulation to 
be Mach ^ 1.2, which, using the |Ostriker (1999) formula 



into Equation [S] results in a rather fast orbital decay with 
Tod{r = 20 pc) ^ 1 Myr, quite consistent with the orbital 
evolution of the BH seen in Figure |3] When the black hole 
fall within the innermost r = 5 pc of the nuclear disk, its 
relative velocity drop below 100 km/s while the local sound 
speed reach Cs = 310 km/s. In this strongly subsonic regime, 
the drop of /^sas)j^^ ^ ^ ^ q compensated by 

the increase of the ratio p/c^ and the orbital decay timescale 
rises to Tod{r ^ 1 pc) 10 Myr. This, together with the 



fact that the gas mass present between the two black holes 
is comparable or smaller than the sum of the masses of the 
two holes below pc scale separations, explains in principle 
why the decay becomes inefficient. We will return on this 
point in the Discussion below. While the BH is sinking to- 
wards the center of the nuclear disk, its orbit is getting more 
and more circular, as found bylDotti et al.|(2006, 20071. This 



will reduce the relative velocity of the BH with respect to 
the gas disk even more, and strengthens the robustness of 
our conclusion, namely that BH pairing in a stable, pressure- 
supported, nuclear disk is a rather slow mechanism. 

The issue of gap-opening is an important point in the 
process of black hole binary formation and its final coales- 
cence. If the black holes are able to open a circumbinary 
gap, the density of the gas surrounding the black holes drops 
significantly and dynamical friction is no longer efficient to 
make the binary separation decrease even further. Decay 
might continue if there are other viscous processes acting in 
the disk, as in the case of migration in protoplanetary disks, 
but the timescale will be related to the specific nature of 
the viscous process (e.g. transport of angular momentum via 
MHD or density waves supported by the self-gravity of the 
gas) which we either do not include (e.g. MHD instabilities) 
or cannot model properly at scales near our resolution limit 
(e.g. density waves). One can determine the gap-opening cri- 
terion by comparing the gap-closing and gap-opening times 
as it has been done in the context of planetary rings ( ^Gol-| 



dreich Tremaine | [19821 . In our nuclear disk of thickness 



60 pc, the black hole mass gap-opening criterion can be 
written ( Escala et al.||200 5| eq. 6) 



Mb 



Ar/h 



47ryW/(g-)(yW) 



M(< r) 



(9) 



When the two black holes fall into the innermost 5 pc, 
M 0.15 and f^^^'^M) ^ 0.05. As shown in the previ- 
ous section, the black hole should clear a gap of typical size 
Ar ^ 2.5 pc to prevent the formation of a hydrodynami- 
cal wake and make the dynamical friction inefficient. Con- 
sequently, the black hole critical mass is 2 x 10^ M©, well 
above the mass of our black hole particle which could never 
open a gap in such a nuclear disk. Even if a Bondi-Hoyle 
accretion modeled was taken into account in the simulation, 
^Callegari et al.| ( |2011) showed that, in minor merger simula- 
tions, one of the black hole could grow by nearly an order of 
magnitude in mass at the most, which would not be enough 
for our 2.5 X 10^ Mq black hole to reach the critical mass. 
However, since the observed black hole masses range from 
10^ to a few 10^ M©, the formation of a circumbinary gap 
by a more massive black hole is possible in our disk. 



4 CONCLUSIONS AND DISCUSSION 

Our numerical study is the first one that uses an AMR code 
to address the formation of a nuclear gas disk resulting from 
an equal-mass dissipative galaxy merger and the concurrent 
formation of a SMBH binary at parsec scales within such a 
disk. Our results confirm previous findings of SPH simula- 
tions with the GASOLINE code, with particle splitting that 



started from very similar initial conditions ( [Mayer et al. 
2007| [Mayer k Kazantzidis||2008|). After the initial asym- 



metric black hole injection following the final merger of the 
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galactic cores, the dynamical friction from the gaseous back- 
ground dominates the one due to the stellar background and 
makes the black hole separation falls from 100 pc to a few 
parsecs in less than 10 Myr, in qualitatively good agree- 



The exact timescale of the sinking of the SMBH binary 
is however dependent on resolution, as was also noticed in 



ment with [Mayer et al. ( 2007p results. Indeed in the refer- 
ence simulation adopting a polytropic equation of state with 
Fd = 7/5, the mass, size and typical density of the nuclear 
disk that arises (respectively, ^ 10^ and 60 pc) are ex- 
tremely close to those of the nuclear disk in Mayer et al.^ 
(2007), as is its characteristic temperature. Due to the poly- 



tropic equation-of-state we adopted for the high-density gas 
(Fd = 7/5), the nuclear disk is hot, its thickness is ^ 60 pc 
consistently with observed nuclear regions in ULIRGs and 
the flow is laminar (the turbulence generated by the colli- 
sion of the galaxy cores is rapidly damped). Being embed- 
ded in a background with properties similar to that in the 
[Mayer et al.| (|2007|) calculations, in particular with compa- 
rable density and temperature, the pair of SMBHs sinks by 
dynamical friction very rapidly after the merger, and be- 
comes bound also at a similar separation, about 6 pc. We 
have also run two lower resolution simulations with different 
polytropic indexes, one with Fd = 7/5 as in the reference run 
and another one with Fd = 5/3; the comparison shows the 
same trend found by [Mayer et al. ( 2007 ) , namely that with 
a stiffer EOS (Fd = 5/3) the orbital decay is slower, both 
because the decay regime is more subsonic and because the 
gaseous background has a lower characteristic density, both 
effects going in the direction of reducing dynamical friction. 
At a closer inspection, there are however some differ- 



ences between the SPH nuclear disks of Mayer et al. ( 2007 ) 



and that in the new RAMSES simulations presented in this 
paper. First, the disk in the SPH simulations displayed a 
stronger spiral structure since its appearance, as a result 
of a higher central surface density and thus a stronger self- 
gravity (the Toomre parameter was close to 1.5 while ^2.0 
in our disk, see top row of fig. [2]). The stronger spiral modes 
lead to a more effective transfer of angular momentum and 
thus a continuously increasing central density. Second, the 
level of non-thermal motions in the gas, usually termed "tur- 
bulence" , was a factor of about 3 higher than in the RAMSES 
run (^ 100 km/s instead of 30 — 40 km/s), likely as a re- 
sult of the stronger self- gravitating response. The reason for 
these differences are not clear. The higher numerical viscos- 
ity in the SPH runs, due to the use of explicit Monaghan 
artificial viscosity, especially the quadratic term in shock dis- 
sipation during the merger, is expected to dissipate gas mo- 
tions into thermal energy. This would damp the turbulence 
faster, at odds with the higher turbulence seen in the SPH 
simulation, but would also generate viscous transport of an- 
gular momentum, leading to an increase in central density 
and therefore a higher susceptibility to gravitoturbulence, 
at least in the innermost region. Another possible cause of 
the difference is the fact that in the GASOLINE runs of Mayer 



et al. ( 2007| the energy equation was solved, including shock 



dissipation via artificial viscosity, possibly producing a nu- 
clear disk with lower entropy immediately after the merger 
relative to the RAMSES runs adopting a fixed polytropic equa- 
tion of state which does not account, by construction, for 
entropy dissipation in shocks. Finally, as we explain below, 
resolution in SPH and AMR runs is not guaranteed to be 
identical even if the setup is designed to be as close as pos- 
sible as in the case that we are discussing. 



Mayer & Kazantzidis (2008), in particular becomes shorter 



for increasing resolution. At higher resolution the dynamical 
friction wake is better resolved as density gradients are bet- 
ter captured, this being a likely reason behind the faster dy- 
namical friction timescale (in collisionless systems an analo- 
gous resolution dependence of dynamical friction, mediated 
by gravitational softening, has been reported extensively in 
the literature, see e.g. Colpi et al. (1999). The formation 



of a front shock and a trailing hydro dynamical wake which 
exerts a gravitational drag on the black hole, the so called 
dynamical friction wake, is a remarkable result of our RAMSES 
simulations. The RAMSES simulations presented in this paper 
are indeed the first that allow to capture the wake of dy- 
namical friction very clearly in a highly dynamic situation 
such as that of a nuclear disk arising from merging galaxies. 
We show that this gravitational drag is due to overdense gas 
within 2 — 3 pc behind the black hole and that the efficiency 
of the hydrodynamical friction peaks at transonic regime, in 
fair agreement with the analytical prediction from Ostriker] 
(1999). Although we did not prove that the dynamical fric- 
tion wake structure has already converged, the fact that it 
is very clearly resolved, while it was not in previous cal- 
culations conducted with SPH simulations, prompts us to 
believe that our results are quantitatively more robust than 
those in previous works. 

Our RAMSES simulations, owing to the aggressive re- 
finement enabled by the AMR technique, allow to reach a 
spatial resolution of 0.1 pc in the center of the disk where 
the SMBHs are sinking, which is ten times better than the 



nominal resolution in,Mayer et al. (2007). We caution, how- 



ever, that comparing the resolution in AMR and SPH is 
not straightforward. In particular, it in the SPH simulations 
of Mayer et al. (2007) and Mayer Kazantzidis (2008^ we 
adopted a fixed gravitational softening in the high resolution 
region after particle splitting, which can thus be considered 
the actual limit of spatial resolution in those calculations 
(the SPH smoothing length being smaller in high density 
regions due to the large number of particles employed). In 
RAMSES the gravitational force resolution is not fixed, rather 
it is tied to the cell size, therefore it shrinks as the refinement 
is applied to the cells , a situation more reminiscent of what 
happens with SPH codes adopting an adaptive softening 
length (e.g. [Bate k Burkert|1997 ). We notice that in Mayer[ 



& Kazantzidis 



( 2008 ) a run was presented with a maximum 



spatial resolution, in terms of gravitational softening, com- 
parable to the highest resolution simulation presented here; 
although the simulation was carried out only for a few orbits 
after the black holes have formed a binary rather than for 
many orbits as in the RAMSES simulations presented here, in 
both cases the separation of the SMBHs appears to fluctu- 
ate significantly, with no clear signs of sustained decay be- 
low parsec separation. Below such separation it is expected 
that dynamical friction will become inefficient because the 
mass of the gaseous background enclosed within the orbit 
of the two SMBHs becomes smaller than the mass of the 
SMBH binary. In addition, in the innermost 5 pc of the nu- 
clear disk, the sound speed reaches ^ 300 km/s, a subsonic 
regime for the motion of SMBHs, which also implies an inef- 
ficient dynamical friction. In polytropic equilibrium models 



of nuclear disks, Escala et al.| ([2004| reported an asymmet- 
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ric torque due to an ellipsoidal deformation of the density 
distribution around the binary, which extracted angular mo- 
mentum from the binary allowing it to shrink further. This 
is not seen here nor in the previous SPH simulations of some 
of us. One reason might be that the thermodynamical con- 
ditions and density profile of the disk that develops here are 
different from those in equilibrium disk models. In particu- 
lar, the system analyzed by Escala et al. (2004) had a sound 



speed of about 60 km/s even at the center, hence it was a 
much colder gas disk than ours. On the other hand, a colder 
disk could act against the decay in two ways, namely by 
driving a stronger mass inflow by gravitational torques that 
steepens the central density on a dynamical timescale and 
might thus evacuate the gas around the black hole binary, or 
might provide more favourable conditions for the opening of 
a gap by the binary after it has become bound. Finally, it is 
important to note that binary decay stalls around 2 — 3 pc, 
a separation at which the orbital evolution might not be 
correctly captured with only ^ 20 AMR cells. A similar res- 
olution problem, as well as possible associated issues with 
orbit integration accuracy, might have been present also in 
the previous SPH calculations. 

Our findings thus suggest that, as in the case of purely 
stellar backgrounds, a continued decay towards the gravita- 
tional wave regime is not automatically achieved in a gaseous 
background. However, concluding that there is a last parsec 
problem in gaseous backgrounds is highly premature. First 
of all, while the strength of our models, relative to other 
studies adopting idealized nuclear disks, relies in the realis- 
tic disk conditions inherited by the galaxy merger, there are 
still several important simplifications and omissions in the 



physics at play. First, as in Mayer et al. (2007) and Mayer 
& Kazantzidis (2008|), we considered a single phase medium 



described by an effective EOS. In reality the nuclear disk will 
host a complex multi-phase ISM, with possibly a highly in- 
homogeneous density structure (e.g. |Wada k, Norman|2001 ). 
Star formation and supernovae feedback will provide both 
diffuse and localized heating sources, and if the black holes 
are active as AGNs during one or more phases of the merger 
they would likely change the initial conditions of the nu- 
clear disk arising after the merger (i.e. change its density 
and thermal structure, both being crucial for dynamical fric- 
tion). Ongoing work with a new multi-phase ISM, star for- 
mation and feedback scheme, for the moment implemented 
only in GASOLINE , will soon provide a clue on the impor- 
tance of such complexity in the sinking rate of the SMBH 
binary (Roskar et al., in preparation). Furthermore, allow- 
ing for star formation to happen in the nuclear disk, another 
missing ingredient, will have an impact since stars can aid 
the decay in regions where the gas becomes inefficient, as 
long as stars move on sufficiently non circular orbit with re- 
spect to the frame of the binary, in order to keep the loss 



cone continuously filled ( ^Preto Am aro-S eoane^2010( Khan 
et al.|20lT ). Indeed in a massive self-gravitating disk such as 



the one obtained here, stars would likely exhibit centrophilic 
orbits as a result of non-axisymmetric instabilities. One can 
image thus a multi-stage decay, in which gas is the leading 
drag source down to parsec scales, owing to its very effi- 



cient action demonstrated here (Mayer et al. 2007), and stars 



might take over at smaller separations. Future calculations 
capable of capturing both a realistic multi-phase gaseous 
medium and the full stellar dynamical response including 



three-body interactions between the binary and the stars at 
small scales mark the necessary next frontier that will be 
needed to progress further in understanding the shrinking 
of SMBH binaries well below parsec scales. 
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